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Recent  technologies  have  led  to  stringent  control  system  requirements.  This  has  increased  the 
importance  and  complexity  of  the  analysis  and  design  of  control  systems.  For  example,  the  design 
of  such  control  systems  often  requires  the  solution  of  systems  of  nonlinear  equations  of  high  order. 
Homotopy  algorithms  which  solve  systems  of  nonlinear  equations  form  the  basis  of  the  software 
package  known  as  HOMPACK.  However,  because  of  the  high  dimension  and  structure  of  modem 
control  problems,  the  HOMPACK  code  cannot  be  applied  directly  to  many  of  the  problems  which 
arise  in  control  system  design  and  analysis.  Thus,  the  primary  objective  of  the  current  research 
is  to  extend  and  develop  homotopy  algorithms  for  a  variety  of  computational  problems  in  control. 
In  addition,  these  problems  are  being  examined  in  the  context  of  the  algebraic  and  differential 
geometry  on  which  the  homotopy  methods  are  based.  This  will  enable  the  classification  of  solutions 
to  a  particular  problem  and  consequently  will  allow  the  analyst  or  designer  to  extract  the  most 
desirable  solution. 


Since  the  funding  arrived  too  late  to  recruit  a  student  for  Fall  1989,  work  was  begun  on 
preconditioned  conjugate  gradient  algorithms  for  solving  large,  sparse  systems  of  linear  equations. 
Regardless  of  what  control  algorithms  are  eventually  developed,  they  will  depend  on  sparse  matrix 
technology.  Good  progress  has  been  made  in  this  regard,  as  numerous  experiments  on  realistic  large 
problems  have  been  conducted  with  various  PCG  algorithms,  GMRES,  and  Orthomin(k).  This 
work  has  been  submitted  to  the  SIAM  Journal  on  Optimization,  and  current  analyses  of  SYMLQ 
will  be  presented  at  the  Copper  Mountain  Conference  on  Iterative  Methods  in  April,  1989. 

In  January,  1990,  a  graduate  student,  Dragan  2igic,  began  work  on  the  optimal  projection 
equations  in  papers  of  Bernstein,  Hyland,  and  Richter.  After  some  background  study  in  linear 
systems  theory  and  optimal  control,  Zigic  should  be  able  to  make  some  progress  on  the  control 
design  problems  described  in  Section  2  of  the  proposal. 

Fifteen  conference  presentations  were  given  in  1989.  Those  published  in  proceedings  were: 

•  Two  on  parallel  homotopy  algorithms  for  a  hypercube  -  A.  Chakraborty,  D.  C.  S.  Allison,  C. 
J.  Ribbens,  and  L.  T.  Watson,  “Parallel  orthogonal  decompositions  of  rectangular  matrices  for 
curve  tracking  on  a  hypercube”,  in  Proc.  Fourth  Conf.  on  Hypercube  Concurrent  Computers 
and  Applications,  J.  Gustafson  (ed.),  ACM,  Monterey,  CA,  1989;  A.  Chakraborty,  D.  C. 
S.  Allison,  C.  J.  Ribbens,  and  L.  T.  Watson,  “Parallel  unit  tangent  vector  computation  for 
homotopy  curve  tracing  on  a  hypercube”,  in  Proc.  1990  ACM  Eighteenth  Annual  Computer 
Science  Conference,  Washington,  DC,  1990, 103-108. 


•  Two  optimal  control  problems  -  G.  Vasudevan,  L.  T.  Watson,  and  F.  H.  Lutze,  “A  homotopy 
approach  for  solving  constrained  optimization  problems”,  in  Proc.  Amer.  Control  Conf, 
Pittsburgh,  PA,  1989,  780-785;  G.  Vasudevan,  F.  H.  Lutze,  and  L.  T.  Watson,  “A  homotopy 
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method  for  space  flight  rendezvous  problems”,  in  Proc.  AAS/AIAA  Astrodynunics  Specialist 
Conf.,  AAS,  Stowe,  VT,  1989. 

•  Two  chapters  in  books  -  L.  T.  Watson,  “Modem  homotopy  methods  in  optimization”,  in 
Impacts  of  Recent  Computer  Advances  on  Operations  Research,  R.  Sharda,  B.  L.  Golden,  £. 
Wasil,  0.  Bald,  W.  Stewart  (eds.),  North-HoUand,  New  York,  1989,  555-565;  L.  T.  Watson 
and  M.  P.  Kamat,  “Homotopy  algorithms  for  engineering  analysis”,  in  Supercomputing  in 
Engineering  Analysis,  H.  Adeli  (ed.),  Marcel  Dekker,  New  York,  1990. 

Journal  papers  completed  and  submitted  were: 

•  A  fluid  mechanics  application  -  C.  J.  Ribbens,  C.  Y.  Wang,  L.  T.  Watson,  and  K.  A.  Alexander, 
“Vortidty  induced  by  a  moving  elliptic  belt”,  Comput.  &  Fluids. 

•  Two  papers  on  truss  design  via  homotopy  methods  -  V.  Arun,  C.  F.  Reinholtz,  and  L.  T.  Wat¬ 
son,  “Enumeration  and  analysis  of  variable  geometry  tmss  manipulators”,  ASME  J.  Mecha¬ 
nisms,  Transmissions  Automation  Design;  V.  Arun,  C.  F.  Reinholtz,  and  L.  T.  Watson,  “New 
homotopy  solution  techniques  applied  to  variable  geometry  trusses”,  ASME  J.  Mechanisms, 
Transmissions  Automation  Design. 

•  A  paper  on  parallel  curve  tracking  -  A.  Chakraborty,  D.  C.  S.  Allison,  C.  J.  Ribbens,  and  L. 
T.  Watson,  “Unit  tangent  vector  computation  for  homotopy  curve  tracking  on  a  hypercube”. 
Parallel  Comput. 

•  A  preliminary  study  of  linear  algebra  techniques  applicable  to  large  sparse  control  problems  - 
K.  M.  Irani,  C.  J.  Ribbens,  H.  F.  Walker,  L.  T.  Watson,  and  M.  P.  Kamat,  “Preconditioned 
conjugate  gradient  algorithms  for  homotopy  curve  tracking”,  SIAM  J.  Optim. 

•  A  solid  mechanics  application  -  C.  Y.  Wang  and  L.  T.  Watson,  “Rotation  of  polygonal  space 
structures”,  J.  Astronaut.  Sci. 

Period:  3/1/90  -  2/28/91 

During  the  spring  and  summer  Dragan  ^gic  worked  through  most  of  Kwakemaak  and  Sivan’s 
optimal  control  book,  read  several  homotopy  papers,  and  studied  the  optimal  projection  papers  of 
Bernstein,  Collins,  Hyland,  and  Richter  in  depth.  The  collaboration  with  Bernstein  and  Collins 
of  Harris  Corporation  in  Melbourne  has  been  extensive  and  very  productive.  They  have  provided 
us  with  test  data  and  guidance  as  to  what  is  important,  and  we  have  improved  their  theoretical 
results  and  numerical  algorithms.  The  first  fruit  of  this  collaboration  is  a  conference  paper  entitled 
“A  homotopy  method  for  solving  the  optimal  projection  equations  for  the  reduced  order  model 
problem”  to  be  presented  at  the  IEEE  Sontheastcon  meeting  in  Williamsburg  in  April  1991.  This 
work,  also  part  of  Zinc’s  M.S.  thesis,  is  summarized  here: 

The  optimal  projection  approach  is  utilized  on  various  problems  arising  in  optimal 
control.  Hyland  and  Bernstein  [1]  ^ve  theoretical  results  for  the  application  of  that 
method  to  the  reduced  order  model  problem,  which  is  to  find  a  reduced  order  model 

®m(t)  —  Am  ®m(0  -Sm  **(0> 

Vmif)  ~  Cfn^mit), 


for  the  system 


x{t)  =  Ax{t)  +  Bri{t), 
y{t)  =  C  x{t). 
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that  minimizes  the  quadratic  model  reduction  error 


=  ^hm  ^  (^(0  “  ym(0)]  • 

It  is  assumed  that  both  the  system  and  the  model  are  asymptotically  stable,  control¬ 
lable  and  observable.  Necessary  conditions  for  the  optimal  reduced  order  model  can  be 
expressed  in  the  form: 


0  =  t[AQ  +  QA*  +  BVB% 

(1) 

0  =  [A^P  +  PA  +  C^  RC]t, 

(2) 

rank  {Q)  =  rank  (P)  =  rank  (Q  P)  =  Um, 

(3) 

where  Um  is  the  degree  of  the  model,  Q  and  P  are  pseudo-Gramians  (analogous  to  Grami- 
ans  and  have  rank  deficiency)  and  the  skew  projection  operator  r  is  a  nonlinear  function 
of  Q  and  P.  The  optimal  model  {Am^BmyCm)  can  be  computed  as  a  nonlinear  function 
of  (A,  B,  C)  and  Q  and  P. 

Equations  (1)  and  (2),  called  modified  Lyapunov  equations,  resemble  standard  matrix 
Lyapunov  equations,  but  are  highly  nonlinear  since  they  contain  r.  The  algorithm  pro¬ 
posed  in  this  paper  utilizes  probability-one  homotopy  theory  as  the  main  tool  for  solving 
the  system  (l)-(3).  There  is  a  family  of  systems  (a  homotopy)  which  make  a  continu¬ 
ous  transformation  from  some  initial  system  to  the  final  system  (l)-(3).  Each  system 
along  the  homotopy  path  is  itself  solved  by  a  homotopy  algorithm-a  homotopy  within 
a  homotopy,  so  to  speak.  The  central  theorem  of  the  paper  shows  the  validity  of  the 
whole  process,  i.e.,  determines  the  class  of  initial  systems  which  certainly  lead  to  the  final 
system  along  a  homotopy  path.  Another,  significantly  simpler,  homotopy  is  used  to  solve 
the  initial  problem.  Finally,  it  is  shown  how  the  optimal  solution  to  the  reduced  order 
model  problem  can  be  computed  in  an  easy  way  from  a  solution  to  the  system  (l)-(3). 

[1]  D.  C.  Hyland  and  D.  S.  Bernstein,  The  Optimal  Projection  Equations  for  Model 
Reduction  and  the  Relationships  Among  the  Methods  of  Wilson,  Skelton  and  Moore, 
IEEE  Transactions  on  Automatic  Control,  Vol.  AC-30,  No.  12,  December  1985,  pp. 
1201-1211. 

Nine  conference  presentations  were  given  in  1990.  Those  published  in  proceedings  were: 

•  Two  on  parallel  homotopy  algorithms  for  a  hypercube — A.  Chedcraborty,  0.  C.  S.  Allison,  C. 
J.  Ribbens,  and  L.  T.  Watson,  “Low  dimensional  homotopy  curve  tracking  on  a  hypercube”, 
in  Ptoc.  1990  Internat.  Cord,  on  Parallel  Processing,  Vol.  HI,  P.-C.  Yew  (ed.),  St.  Charles, 
IL,  1990,  44-51;  A.  Chakraborty,  D.  C.  S.  Allison,  C.  J.  Ribbens,  and  L.  T.  Watson,  “Parallel 
homotopy  curve  tracking  on  a  hypercube”,  in  Parallel  Processing  for  Scientific  Compating,  J. 
Dongarra,  P.  Messina,  D.  C.  Sorensen,  and  R.  G.  Voigt  (eds.),  SIAM,  Philadelphia,  PA,  1990, 
149-153. 

•  Two  mechanisms  problems — V.  Arun,  C.  F.  Reinholtz,  and  L.  T.  Watson,  “Application  of 
new  homotopy  continuation  techniques  to  variable  geometry  trusses”,  in  Cams,  Gears,  Robot 
and  Mechanism  Design,  DE-Vol.  26,  A.  Pisano,  M.  McCarthy,  S.  Derby  (eds.),  ASME,  New 
York,  1990,  87-92;  V.  Arun,  C.  F.  Reinholtz,  and  L.  T.  Watson,  “Enumeration  and  analysis 
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of  variable  geometry  truss  manipulators”,  in  Cams,  Gears,  Robot  a,ad  Mecbaaism  Design, 
DE-Vol.  26,  A.  Pisano,  M.  McCarthy,  S.  Derby  (eds.),  ASME,  New  York,  1990,  93-98. 

•  A  NATO  ASI  invited  lecture — L.  T.  Watson,  “Numerical  analysis  of  nonlinear  equations  in 
computer  vision  and  robotics”,  in  Numerical  Linear  Algebra,  Digital  Signal  Processing  and 
Parallel  Algorithms,  NATO  ASI  Series  F,  G.  Golub  and  P.  Van  Dooren  (eds.),  Springer- Verlag, 
Berlin,  1990,  695-704. 

Journal  papers  completed  and  submitted  were: 

•  Four  solid  mechanics  applications — L.  T.  Watson  and  C.  Y.  Wang,  “Large  deformations  of 
rotating  polygonal  space  structures”,  Comput.  Math.  AppL;  C.  Y.  Wang  amd  L.  T.  Watson, 
“Large  deformations  of  a  whirling  elastic  cable”,  Acta  Mecb.;  J.  Rakowska,  R.  T.  Haftka,  and 
L.  T.  Watson,  “An  active  set  algorithm  for  tracing  parametrized  optima”,  Struct.  Optim.;  J. 
Rakowska,  R.  T.  Haftka,  and  L.  T.  Watson,  “Multi-objective  control-structure  optimization 
via  homotopy  methods”,  SIAM  J.  Optim. 

•  A  fluid  mechanics  problem  in  biology — Zs.  Nagy-Ungvarai,  J.  J.  Tyson,  S.  C.  Muller,  L. 
T.  Watson,  and  B.  Hess,  “Experimental  study  of  spiral  waves  in  the  Ce-catalyzed  Belousov- 
Zhabotinskii  reaction” ,  J.  Pbys.  Cbem. 

•  Three  survey  papers  on  optimization — L.  T.  Watson  and  A.  P.  Morgan,  “Serial  and  parallel 
global  optimization  of  polynomial  programs  via  homotopy  algorithms”,  SIAM  J.  Optim.;  L. 
T.  Watson,  “Globally  convergent  homotopy  algorithms  for  nonlinear  systems  of  equations”. 
Nonlinear  Dynamics;  L.  T.  Watson,  “A  survey  of  probability-one  homotopy  methods  for  en¬ 
gineering  optimization”,  Arabian  J.  Sci.  Engrg. 

•  A  paper  on  sparse  matrix  technology — C.  deSa,  K.  M.  Irani,  C.  J.  Ribbens,  L.  T.  Watson, 
and  H.  F.  Walker,  “Preconditioned  iterative  methods  for  homotopy  curve  tracking”,  SIAM  J. 
Sci.  Stat.  Comput. 

•  An  analysis  of  parallel  homotopy  algorithms — A.  ChaJiraborty,  D.  C.  S.  Allison,  C.  J.  Ribbens, 
and  L.  T.  Watson,  “Analysis  of  function  component  complexity  for  hypercube  homotopy 
algorithms”,  IEEE  Trans.  Parallel  Distrib.  Sys. 

•  A  theoretical  numerical  analysis  paper — G.  Ulrich  and  L.  T.  Watson,  “Positivity  conditions 
for  quartic  polynomials”,  SIAM  J.  Sci.  Stat.  Comput. 

•  Another  survey  paper  representing  a  major  application  of  homotopy  methods  to  circuit 
simulation — R.  C.  Melville,  Lj.  Trajkovic,  S.-C.  Fang,  and  L.  T.  Watson,  “Globally  convergent 
homotopy  methods  for  the  DC  operating  point  problem”,  SIAM  J.  Optim. 

HOMPACK  has  become  the  basis  for  a  major  circuit  simulation  code  development  eflort  within 
AT&T  Bell  Labs.  Over  200  requests  for  HOMPACK  have  been  received,  indicating  ,nat  the  work 
is  having  an  impact,  thanks  to  AFOSR  support. 

Period:  3/1/91  -4/30/92 

The  “homotopy  within  a  homotopy”  scheme  described  in  the  previous  section,  similar  to 
the  homotopy  algorithms  of  Collins  and  Richter,  turned  out  to  have  theoretical  and  computationad 
flaws.  Zigic  devised  a  smooth  homotopy  based  on  the  Drazin  inverse  {QP)*  of  i^P,  and  successfully 
solved  a  number  of  problems  with  this  homotopy.  This  Drazin  inverse  based  homotopy  is  described 
in  the  Williamsburg  IEEE  Southeastcon  proceedings  (April,  1991),  and  in  a  chapter  in  Zigic’s  MS 
thesis. 
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In  late  spring  of  1991  Zigic  discovered  a  new  formulation  of  (l)-(3),  based  on  the  contragredient 
transformation,  that  does  not  involve  the  explicit  calculation  of  Q  and  P  or  any  generalized  inverse. 
A  homotopy  based  on  this  contragredient  formulation  successfully  solved  every  problem  we  could 
find  in  the  literature,  as  well  as  a  realistic  space  structure  control  problem  from  Marshall  Space 
Flight  Center  in  Huntsville.  This  algorithm  became  the  basis  for  three  journal  papers  (listed  below) 
and  Zigic’s  MS  thesis. 

Zigic  spent  the  summer  of  1991  working  on  technical  details  for  choosing  a  good  starting  point, 
and  submitted  this  work  for  the  First  Conference  on  Control  Applications  to  be  held  in  Dayton  in 
September  1992.  This  paper  was  accepted,  and  ranked  among  the  top  10%  of  all  papers  submitted. 
Zigic  graduated  and  returned  to  Yugoslavia  in  August  1991. 

Yuzhen  Ge,  a  Ph.D.  in  mathematics,  took  over  from  Zigic  in  fall  1991.  She  spent  six  months 
reading  Zigic’s  thesis,  homotopy  papers,  and  reports  by  Bernstein,  Collins,  Hyland,  and  Richter. 
WhUe  the  Zigic  algorithm  is  accurate  and  robust,  it  involves  too  many  unknowns  to  be  practical 
for  large  scale  problems.  Ge  has  spent  all  of  1992  implementing  and  testing  an  input  normal  form 
homotopy  suggested  by  Collins,  which  has  a  very  small  number  of  unknowns  and  might  be  practical 
for  large  scale  problems.  That  work  is  reported  in  the  attached  document  “An  input  normal  form 
homotopy  for  the  optimal  model  order  reduction  problem”.  Unfortunately  the  input  normal 
form  equations  are  inherently  unstable,  so  a  practical,  robust  homotopy  for  large  scale  problems 
remains  an  open  question. 

Since  the  beginning  of  this  grant,  two  MS  students  and  one  Ph.D.  student  have  been  supported. 
The  MS  theses  are 

•  Kashmira  M.  Irani,  “Preconditioned  sequential  and  parallel  conjugate  gradient  algorithms  for 
homotopy  curve  tracking,”  M.S.  thesis.  Dept,  of  Computer  Sci.,  Virginia  Polytechnic  Institute 
and  State  Univ.,  Blacksburg,  VA,  May  1990. 

•  Dragan  Zigic,  “Homotopy  methods  for  solving  the  optimal  projection  equations  for  the  reduced 
order  model  problem,”  M.S.  thesis.  Dept,  of  Computer  Sci.,  Virginia  Polytechnic  Institute 
ajid  State  Univ.,  Blacksburg,  VA,  June  1991. 

Fourteen  conference  presentations  were  given  during  the  current  period  (3/1/91-4/30/92): 

•  Fifth  SIAM  Conference  on  Parallel  Processing  for  Scientific  Computing,  Houston,  TX,  March, 
1991  (2  papers). 

•  IEEE  Southeastcon,  Williamsburg,  VA,  April,  1991. 

•  Sixth  Distributed  Memory  Computing  Conference,  Portland,  OR,  April,  1991. 

•  Second  International  Conference  on  Industrial  and  Applied  Mathematics,  Washington,  DC, 
July,  1991  (4  papers). 

•  Computational  Structures  Technology,  Edinburgh,  Scotland,  August,  1991  (keynote  lecture). 

•  1991  International  Conference  on  Parallel  Processing,  St.  Charles,  IL,  August,  1991. 

•  Fourth  SIAM  Conference  on  Applied  Linear  Algebra,  Minneapolis,  MN,  Sept.,  1991. 

•  Sixth  IIMAS-UNAM  Workshop  on  Numerical  Analysis  and  Optimization,  Oaxaca,  Mexico, 
January,  1992. 

•  Copper  Mountain  Conference  on  Iterative  Methods,  Copper  Mountain,  CO,  April,  1992. 

•  Scalable  High  Performance  Computing  Conference,  Williaunsburg,  VA,  April,  1992. 
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Those  conference  presentations  published  in  refereed  proceedings  were: 

•  L.  T.  Watson,  “Numerical  analysis  of  nonlinear  equations  in  computer  vision  and  robotics”, 
in  Numericai  linear  Algebra,  Digital  Signal  Processing  and  Parallel  Algorithms,  NATO  ASI 
Series  F,  Computer  and  Systems  Sciences,  Vol.  70,  G.  Golub  and  P.  Van  Dooren  (eds.). 
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AN  INPUT  NORMAL  FORM  HOMOTOPY  FOR  THE  OPTIMAL 
MODEL  ORDER  REDUCTION  PROBLEM 

Yuzhen  Ge^  Emmanuel  G.  Collins,  Jr.*,  Layne  T.  Watson^ 

1.  Introduction. 

The  optimal  model  reduction  problem,  i.e.,  the  problem  of  approximating  a  higher  order 
dynamical  system  by  a  lower  order  one  so  that  a  model  reduction  criterion  is  minimized,  is  of 
significant  importance  and  is  under  intense  study.  Several  earlier  attempts  to  apply  homotopy 
methods  to  the  optimal  model  order  reduction  problem  v’«»re  not  entirely  satisfactory,  Richter 
[1],  [2]  devised  a  homotopy  approach  which  only  estimated  certain  crucial  partial  derivatives  and 
employed  relatively  crude  curve  tracking  techniques.  Zigic,  Bernstein,  Collins,  Richter,  and  Watson 
[3],  [4],  [5],  [6]  formulated  the  problem  so  that  numerical  linear  algebra  techniques  could  be  used 
to  explicitly  calculate  partial  derivatives,  and  employed  sophisticated  homotopy  curve  tracking 
algorithms,  but  the  number  of  variables  made  large  problems  intractable.  We  propose  here  several 
ways  to  reduce  the  dimension  of  the  homotopy  map  so  that  large  problems  are  computationally 
feasible. 

The  problem  can  be  formulated  as:  given  the  asymptotically  stable,  controllable,  and  observable 
time  invariant  continuous  time  system 

i(t)  =  Ai(t)  +  Bu(t), 
y{t)  =  C  x(f), 

where  A  €  R’*’'",  B  €  R"^"*,  C  €  the  goal  is  to  find  a  reduced  order  model 

®Tn(f)  —  Am  "b  ■®m  j.2\ 

ym(0  =  G'm  ®m(0> 

where  Am  6  R'*”*’'""*,  Bm  G  R""*'”*,  Cm  €  R*’*”"*,  <  n  which  minimizes  the  cost  function 

=  ^lu^  E  [(j/ ~  J/m)  E(y  — ym)]j  (3) 

where  the  input  u{t)  is  white  noise  with  synunetric  amd  positive  definite  intensity  V  and  R  is  a 
symmetric  and  positive  definite  weighting  matrix. 

The  optimal  projection  equations  of  Hyland  and  Bernstein  [10],  [11],  described  in  Section  5,  are 
basis  independent  and  correspond  to  the  maximum  number  of  degrees  of  freedom  one  could  plausibly 
use.  Richter  and  Collins  [3]  use  this  majdmum  number,  and  Zigic  [4]  reduced  it  somewhat.  At  the 
other  extreme,  the  minimal  number  of  unknowns  corresponds  to  the  input  normal  form  described  in 
Section  2,  and  developed  into  a  homotopy  algorithm  in  Sections  3  and  4.  Subtle  differences  between 
the  optimal  projection  equations  and  input  normal  form  formulations  are  explored  in  Section  5. 
Section  6  gives  numerical  results  for  the  input  normal  form  homotopy  on  the  test  set  of  Zigic  [4]. 


f  Department  of  Computer  Science,  Virginia  Polytechnic  Institute  and  State  University,  Blacks¬ 
burg,  VA  24061-0106. 

*  Harris  Corporation,  Government  Aerospace  Systems  Division,  Melbourne,  FL  32902. 
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2.  Input  normal  form  formulations. 

The  following  theorem  is  n^ded  to  present  the  homotopy  method  for  the  input  normal  form. 
Theorem  1  [7].  Suppose  Am  is  asymptotically  stable.  Then  for  every  minimal  {Am,  Bm,  Cm), 
i.e.,  {Am,Bm)  is  controllable  and  {Am, Cm)  is  observable,  there  exist  a  similarity  transformation 
U  and  a  positive  definite  matrix  ft  =  diag(c>;i,-  •  -  such  that  Am  =  U~^AmU,  Bm  =  Bm, 

and  Cm  =  CmU  satisfy 

0  =  Am  +  Al-\-BmVBl, 


0  =  Aln  +  aAm  +  clRCm. 


In  addition, 


'  (5) 

_  (ClRC„),.-u,(B„VBl)..  _  ^ 

{Am},. - - ,  ,f  5^  uy 

Definition  l.  The  triple  (A„,,5„,Cm)  satisfying  (4)  or  (5)  is  said  to  be  in  input  normal 
form. 

Under  the  assumption  that  a  solution  {Am,  Bm,Cm)  iu  input  normal  form  is  sought,  the  only 
independent  variables  are  Bm  and  Cm,  and  in  this  case  the  domain  is 

{{Am,  Bm,Cm)  •  Am  is  Stable,  {Am,Bm,Cm)  is  minimal  and  in  input  normal  form}. 

The  cost  function  (3)  can  be  written  as 

J{Am,  Bm,Cm)  —  tr  {QR)  (6) 

where  Q  is  a  symmetric  and  positive  definite  matrix  satisfying 

AQ  +  QA^  +  V  =  0,  (7) 


M  R=(  -C^RCm^ 

\0  Am)'  ^  \-ClRC  ClRCm)' 


Q  can  be  written  as 


t)’ 


'  BVB'^  BVBl 
BmVB^  BmVBl 


where  Qx  e  R"’'",  Qx2  e  R'**'”’",  and  Qj  €  R"- 

The  goal  of  minimizing  (6)  under  the  constraints  (4)  and  (7)  leads  to  the  Lagrangian 

B{Am,  Bm,Cm,^,Q)  =  tr[Qfl  +  (A„,  +  Am  +  BmVBm)Mc 

+  {Ain  +  ilAm  +  CiRCm)M„  +  {AQ  +  QA'^  +  7)P] , 

where  the  symmetric  matrices  Mo,  Me,  and  P  are  Lagrange  multipliers. 
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Setting  dLldQ  =  0  gives 


(9) 


3.  A  homotopy  approach. 

A  homotopy  approach  following  [8]  is  now  described.  Let  A/,  Bf,Cf  ,  Rj,  and  Vf  denote  A, 
B,  C,  R,  and  V  in  the  above  and  define 


A(A)  =  Ao  +  A(A/  -  Ao), 
B{X)  =  Bo  +  X{Bf  —  Bo), 
C(A)  =  Co  +  XiCf  -  Co), 


i?(A)  =  Ro  +  X(^Rj  —  Ro), 
V{X)  =  Vo  +  XiVf-Vo). 


For  brevity,  A(A),  5(A),  C(A),  F(A),  zmd  iZ(A)  will  be  denoted  by  A,  B,  C,  V,  and  R  respectively 
in  the  following.  Let 


EbA^^)  =  ^  =  +  P2BJ1V  +  2McB,„V, 

EcA^,X)  =^  =  2R{CmQ2  -  CQn)  +  2RC„,Mo, 


(14) 
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where 


^/Vec  iB^)\ 

-  VVec(C„)y 


VVec(C„)y 

denotes  the  independent  variables  Bm  and  Cm,  Mo  and  satisfy  (12),  and  Q  and  P  satisfy 
respectively  (7)  and  (9).  Vec(P)  for  a  matrix  P  6  is  the  concatenation  of  its  columns: 


Vec(P)  = 


The  homotopy  map  is  defined  as 


and  its  Jacobian  matrix  is 


Define 


^  VVec  [5c.(^,A)]j’ 


Dp(0,X)={Depie,X),Dxpi0,X)). 


HB„{pC\M(J^)=2{P^^^^^B  +  Pi^^Bm)V  +  2Mi^^BmV, 

ScAQ^^\M<<J'>)  =  2R{CmQ^^^  -  CQ'^S)  + 

where  the  superscript  (j)  means  djd0j:  =  §f-.  Using  the  above  definitions,  we  have  for 

=  Bb„  Mp))  +  2(P2  + 

Cf\Bm)kl  /IQN 


and  for  0j  =  {Cm)f 


0{Bm}kl 


=  Jc.  MW)  +  {Q2  +  Mo), 

where  is  a  matrix  of  the  appropriate  dimension  whose  only  nonzero  element  is  e*/  =  1- 
and  can  be  obtained  by  solving  the  Lyapunov  equations 

0  =  +  Ag(^>  + 

0  =  +  P<-»Ji  +  P/4(^>  +  ^ 

Similarly  for  A,  using  a  dot  to  denote  d/dX, 

=  Hb„  {P,  Mo)  +  2Pf2 {BV  +  BV)  +  2(P2  +  M,)P„V, 

=  Pc„ (4  Mo)  +  2RCm (gj  +  Mo)  -2{RC  +  RC)Qn, 


(21) 


where  P  and  Q  are  obtained  by  solving  the  Lyapunov  equations 

o  =  Aq  +  a^  +  +  QAt  +  V, 

0  =  A'^P  +  A^P  +  +  Pii  + 


4.  Numerical  algorithm  for  input  normal  form  homotopy. 

The  initial  point  (^,  A)  =  {Or  0)  =  ((5m)o>(C'm)o»0)  is  chosen  so  that  the  triple 
((•4m)o»(Pm)o5(t7m)o)  is  in  input  normal  form  and  satisfies  p(Oo,Q)  =  0. 

Theorem  3  [9].  Suppose  A  is  asymptotically  stable.  Then  for  every  minimal  {A,B,C),  i.e., 
{A,B)  is  controllable  and  {A,C)  is  observable,  thae  exist  a  similarity  transformation  T  and  a 
positive  definite  matrix  A  =  diag  {di,d2,-  ■  ■,dn)  such  that  A  =  T~^AT,  B  =  T~^B,  and  C  =  CT 
satisfy 

0  =  AA  +  +  BVB'^, 

0  =  A^A  + AA  +  C^PC. 

Definition  2.  The  triple  {A,B,C)  in  the  above  theorem  is  balanced. 

According  to  Moore  [9],  under  certain  conditions,  the  leading  principal  X  block  of  A, 
the  leading  principal  n„,  x  m  block  of  B,  and  the  leading  principal  I  x  rim  block  of  C  in  balanced 
form  are  good  approximations  to  the  reduced  order  model.  This  suggests  that  the  initial  point 
(do,0)  be  chosen  as  follows: 

1)  Transform  the  given  triple  {Af,Bf,Cf)  to  balanced  form  (A(„  Pi,  Cj,). 

2)  Partition  (Ai,  Pi,  Ci)  as 

n 

3)  (Ao,Po»C'o)  is  chosen  as 

4)  The  initial  point  for  the  reduced  order  model  is  chosen  as 

.  _  fyeciBm)o\  _  fyecBi\ 

‘'"■'VVec  {Cm)o  )  ~  {Vec  Cj  ' 

and  (A„,)o  =  An  by  construction. 

5)  Transform  the  initial  point  ((A„,)o,(Pm)o,(^m)o)  to  input  normal  form  so  that  the  initial 

reduced  order  model  is 

Mm)o,(Pm)o,(C,„)o)  =  (r~'(A„.)oT,  T"' (P„)o,  (Cm)oT). 

The  initizJ  point  for  the  homotopy  map  is  then  (^o»0),  where 


Vec  (Pm)o 
Vec  (Cmjo 


)■ 
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Once  the  initial  point  is  chosen,  the  rest  of  the  computation  is  as  follows; 

1)  Set  A  :=  0,  «  := 

2)  Calculate  Am  from  (5),  V,  and  compute  Q  and  P  according  to  (7)  and  (9). 

3)  Evaluate  5  from  (13)  and  Mo  and  Me  according  to  (12). 

4)  Evaluate  the  homotopy  map  p(0,X)  in  (15)  and  Dp{9,X)  in  (16). 

5)  Predict  the  next  point  =  (tf(°^,A(°^)  on  the  curve  7. 

6)  For  k  ;=  0, 1, 2,  •  •  •  until  convergence  do 

Z(k+i)  ^  [Dp(2(‘))]^p(Z(*)), 

where  [JDp(Z)]^  is  the  Moore- Penrose  inverse  of  Dp{Z).  Let  (^i,Ai)  =  lim  Z^'‘K 

k—^oo 

7)  If  Ai  <  1,  then  set  6  :=  A  ;=  Ai,  and  go  to  step  2). 

8)  If  Ai  >  1,  compute  the  solution  0  at  A  =  1.  Am  is  then  obtained  from  (5). 


5.  Compsirison  with  optimal  projection  equations  approach. 

Theorem  4  [10]  [11].  Suppose  {Am,Bm,Cm)  solves  the  problem  (l)-(3).  Then  there  exist 
pseudogramians  Q,  P  that  are  a  solution  to  modified  Lyapunov  equations 

0  =  t[AQ  +  QA*  +  BVB% 

7  :  (22) 

Q  =  [A‘P  +  PA  +  C'RC]t, 


and  satisfy  rank  conditions 


rank  (Q)  =  rank  (P)  =  rank  (Q  P)  =  Um, 
such  that  the  optimal  model  is  given  by 


Am  =  TAG\ 

Bm  —  P  B, 
Cm=CG\ 

where  G  and  P  come  from  a  {G,  M,  P)-factorization  of  QP: 

QP  =  G*MT, 

PG*  =  /n., 


(23) 


(24) 


G,  P  6  R'*’"’*",  M  G  is  positive  semisimple  and  r  =  G*T. 

Equations  (22)  are  called  the  optimal  projection  equations,  which  after  a  lot  of  algebra  described 
in  [5],  can  be  written  in  a  form  suitable  for  computation  as 

UiAWiJiWl  +  EWlA^  +  UiBVB*  =  0,  (n„,  n) 

A^Ul^  +  UlllUiAWi  +  C^  RCWi  =  0,  {num)  (25) 

u,w,-i  =  a.  (»=„) 

The  unknowns  are  Wi  6  R’*’*"”,  Ui  €  R""**”  and  symmetric  E  G 
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Hyland  and  Bernstein  [11]  stated  that  the  optimal  projection  equations  can  have  at  most  |  ^ 

solutions.  It  is  shown  by  the  following  2-dimensional  example  that  this  is  not  true  in  general. 

The  system  [12]  is  given  by 

'‘=(:2;99  -'5200.9) ■ 

Proposition;  For  the  system  (1)  defined  by  (26),  the  solution  set  of  the  optimal  projection 
equations  contains  three  isolated  solutions  and  a  one-dimensional  manifold  parametrized  by  one 
element  of  either  Wi  or  U\. 

Proof.  The  three  isolated  solutions  are 

Am  =  (-0.005004234),  =  (1.0002127),  Cm  =  (1.0002127), 

Am  =  (-4998.0786133),  Bm  =  (100.0001907),  Cm  =  (100.0001907), 

Am  =  (-0.4659163),  Bm  =  (-1.9404824),  Cm  =  (-1.9404824), 

which  were  obtained  by  both  POLSYS  from  HOMPACK  [17]  and  by  a  homotopy  approach  [4]-[6]. 
The  one- dimensional  manifold  of  solutions  corresponds  to 

=  (-0.4851515),  Hm  =  (0.0),  (7,^  =  (0.0), 

which  can  be  derived  directly  from  the  optimal  projection  equations  as  follows. 

Let  Wi  =  Cfi  =  (x3,X4),  and  S  =  ij.  The  optimal  projection  equations  (25)  for  this 

problem  can  be  written  as 

0  =  aiixjxsXs  +  012X1X22:3*5  +  <121*1*4*5  +  022*1*2*4*5 

-I-  O11X1X5  -h  012*2*5  +  (BVB*)uX3  -t-  (BVB*)2iX4> 

2  2 
0  =  011*1X2*3X5  -|-  012*2*3*5  +  021*1*2*4*5  +  022*2*4*5 

-h  021*1*5  +  022*2*5  +  H*)i2*3  +  .5*)22*4> 

0  =  011*1*3*5  012*2*3*5  +  021*1*3*4*5  +  022*2*3*4*5  (27) 

-h  Oii*3*5  +  021*4*5  +  (C*JiC)uXl  (C*iiC')i2*2> 

0  =  011*1*3*4*5  +  012*2*3*4*5  +  021*1*4*5  +  022*2*4*5 

+  Oi2*3*5  +  022*4*5  +  (C*i2(7)21*l  +  (C'*2iC)22*2> 

0  =  *1*3  -I- *2*4  -  1- 

The  triple  (Am,  Bm,  Cm)  is  given  by 


=  *1(011*3  +  021*4)  +  *2(012*3  +  022*4), 

Bm  =  (*3  *4)  ^  ~  ^>11*3  +  i»21*4, 


Cm  =  CG^  =  (cii  C12) 


ft)- 


Cll*l  +  Cl2*2, 


(28) 


where  T  =  Ui  and  G  =  W^.  Substituting  (26)  into  (27)  and  (28),  setting  Bm  =  is  +  IOO14  =  0  and 
Cm  =  ii  +  100x2  =  0  gives  Xi  =  —  IOO12,  X3  =  —  IOOX4,  and  Am  =  —4852x2X4.  Equations  (27) 
become 


0  =  485200x2X4X5  -  0.49x2x5, 

(29) 

0  =  485200x2X4X5  -  0.49x4X5, 

(30) 

0  =  4852X2X4X5  +  4901x2x5, 

(31) 

0  =  4852x2X4X5  -b  4901x4X5, 

(32) 

0  =  10001x2x4  -  1. 

(33) 

If  X2  =  0  or  X4  =  0,  equation  (33)  wiU  not  be  satisfied.  Only  the  situation  that  X2  ^  0  and  147^0 
is  possible.  Then  equations  (29)-(33)  can  be  reduced  to 

0  =  485200x2X4X5  -  0.49x5, 

0  =  4852x2X4X5  +  4901x5,  (34) 

0  =  10001x2x4  -  1. 


If  xs  ^  0  then  (34)  becomes 

0  =  485200X2X4  -  0.49, 

0  =  4852x2X4  +  4901,  (35) 

0  =  10001x2x4  -  1, 

which  does  not  have  a  solution. 

Thus  X5  =  0,  and  equation  (34)  reduces  to 

10001x2X4  -1  =  0, 


which  gives  Am  =  -4852/10001  =  -0.4851515  corresponding  to  a  one-dimensioncd  manifold 
parametrized  by  X2  or  X4.  Hence  the  solution  Am  =  —0.4851515,  Bm  =  0  and  (7^  =  0  (which  is  not 
controllable  or  observable)  corresponds  to  a  one-dimensional  manifold  of  solutions  of  the  optimal 
projection  equations.  Q.  E.  D. 

The  set  of  solutions  of  the  input  normal  form  equations  contains  the  same  set  of  isolated  solutions 

as  the  optimal  projection  equations,  and  idso  a  fourth  isolated  solution  given  by /Ito  =  Bm  =  Cm  =  0. 

Therefore  the  solution  sets  of  the  two  formulations  are  different. 

The  input  normal  form  equations  can  be  rewritten  as 


0  =  2{P^^B  -f-  P2Bm)V  -1-  2M,BmV, 
0  =  2R{CmQ2  -  CQu)  -I-  2RCmMo. 


Setting  Bm  =  Cm  =  0,  the  equations  become 


0  =  P^iBV, 
0  =  RCQu, 


(36) 


(37) 
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where  Pn  and  Qu  satisfy  respectively 


0  =  Pi2  +  P^Am, 
0  =  AQ 12  +  Ql2-^m, 


which  ha^  a  solution  Pu  =  Qi2 


satisfies 


Am  +  A\^  +  BmV  fim  —  Am  +  =  0 


which  gives  Am  =  0. 


6.  Numerical  results  and  comparisons. 

In  this  section  numerical  results  for  the  input  normal  form  formulations  are  given  for  nine 
systems.  These  systems  have  been  studied  and  solved  in  [4],  [5],  [6]  using  the  optimal  projection 
equations  approach.  Comparison.s  are  made  between  the  input  normal  form  formulations  and  the 
optimal  projection  equations. 

The  cost  J  is  computed  for  each  model  as  tr  {Q  R),  according  to  (38).  For  all  examples 
V  =  iZ  =  /.  AU  the  answers  are  given  in  input  normal  form.  The  solutions  obtained  by  the  input 
normal  form  formulation  are  the  same  as  those  obtained  by  the  optimal  projection  equation  method, 
unless  indicated  otherwise. 

Example  1  [12].  The  system  is  given  by 

lo")- 

The  homotopy  algorithm  converges  to  a  solution  corresponding  to  the  model  of  order  Um  =  1  given 
by 

Am  =  (  -0.00500423 ) ,  Bm  =  { -0.100042 ) ,  Cm  =  {- 10.000021 ) , 

which  is  not  in  the  solution  set  of  [4],  [5],  [6]  by  the  optimal  projection  equation  approach.  This 
model  yields  the  (maximum)  cost  J  =  10000. 

In  the  first  step  of  choosing  an  initial  point,  {A/,Bf,Cf)  is  transformed  to  {Ai„Bb,Cb)i  where 
orthogonal  decompositions  of  two  matrices  are  needed.  If  the  eigenvalues  of  one  of  the  matrices  are 
rearranged  in  ascending  order,  then  a  different  solution  is  obtained,  namely 

Am  =  (  -4998.078625) ,  Bm  =  ( -99.980784) ,  Cm  =  {  -100.019608) . 

This  model  )delds  the  (minimum)  cost  J  =  96.078058. 

Example  2  [13].  The  system  is  given  by 

^=(‘o'  -°io)'  ^=(70  !)' 

A  model  of  order  Um  =  1  is 

An,  =  (-11.979443),  =  ( -4.859135  0.589656),  C,„  =  (2.760762) . 

This  model  yields  the  cost  J  =  0.598377. 
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Example  3  [12].  The  system  is  given  by 

'‘=(-0“  :o"72)'  = 

A  model  of  order  =  1  is 

Am  =  ( -0.838521 ) ,  Bm  =  {  -1.295006 ) ,  Cm  =  { 1.825580 ) . 

This  model  yields  the  cost  J  =  0.107256. 

Example  4  [14].  The  system  is  given  by 

/-I  3  0\  /-2\ 

A=l-1  -1  1  1,  5=1  2j,  C  =  (l  0  0). 

A  model  of  order  =  1  is 

Am  =  ( -0.286334 ) ,  5m  =  ( -0.756748 )  Cm  =  ( 0.878161 ) . 

This  model  yields  the  cost  J  =  1.228834  and  this  solution  is  different  from  that  obtained  by  the 
optimal  projection  equation  method  [4],  [5],  [6].  A  model  of  order  Um  =  2  is 

.  _/ -0.215037  0.753968  ^  ^  _/ 0.655800 \  _/  0.888877  ^ 

“  V  -2.513846  -3.600739  j  ’  "  1,2.683557 )'  “  1,-1.090926 )  ' 

This  model  yields  the  cost  J  =  0.0197781. 

Example  5  [12].  The  system  is  given  by 

/-lO  1  0\  /0\ 

A=  -5  0  1  ,  5=1,  C  =  (l  0  0). 

\-i  0  oy  \lj 

A  model  of  order  nm  =  1  is 


Am  =  (  -0.157898 ) ,  5m  =  (  0.561956 ) ,  Cm  =  (  0.318537 ) . 
This  model  yields  the  cost  J  =  0.0107792.  A  model  of  order  Um  =  2  is 


=(;2: 


-0.139652  0.100607 

-0.600971  -0.448192 


\  R  -  y0.528492\  _  /  0. 

J,  J  ^ 


320438 

.0961019 


This  model  yields  the  cost  J  =  0.000329024. 
Example  6.  The  system  is  given  by 


0  1 

0 

0  \ 

0 

-2  -0.02 

1 

0.01  ' 

0  -  1 

0 

0  0 

0 

1  i 

’  ^  “  0 

0 

0.1  0.001 

-0.1 

-0.001/ 

\0 

1 

C  =  (0  1  0  0). 
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For  this  system,  the  initial u;’s  are  approximately  the  same,  which  leads  to  a  signihcant  numerical 
error  in  computing  Mo  and  the  numerical  failure  of  the  homotopy  algorithm.  Therefore  this  technique 
for  choosing  initial  points  fails,  and  some  modification  to  the  algorithm  is  needed  to  avoid  this  kind 
of  ill  conditioning.  However,  it  is  not  at  all  clear  how  to  systematically  avoid  nearly  equal  cj’s,  and 
this  remains  an  open  question. 

Example  7  [9],  [15].  The  system  is  given  by 


0 

0 

-150 

1 

0 

0 

-245 

0 

1 

0 

-1113 

\o 

0 

1 

-19 

^  =  (0  0  0  1). 


A  model  of  order  rim  =  1  is 


A„  =  (-0.495187),  =  (0.995175),  Cm  =  (0.0148426). 


This  model  yields  the  cost  J  =  4.90749  -10  * .  A  model  of  order  Um  =  2  is 

.  / -0.437964  -0.482612 „  (  0.935911  \  nnnR«onQ7^ 

=  (  2.840074  -3.172419  j  ’  (  -2.518896  j  '  0.00682097 ) . 


This  model  yields  the  cost  J  =  4.159  *10  '.A  model  of  order  rim  =  3  is 

/ -0.437810  -0.483078  -0.0370108  \ 

Am  =  2.826317  -3.135361  -0.612598  , 

\ -4.651841  13.160394  -12.554152/ 


/  0.935746  \ 

.Bm  =  -2.504141  ,  Cm  =  (0.0149143  0.00682180  0.000635413). 

\  5.010819  / 

This  model  yields  the  cost  J  =  4.59  •  10“^°. 

Example  8  [9].  The  system  is  given  by 


(  ° 

1 

0 

0  ^ 

0 

0 

1 

0 

0 

1 

,  B  = 

2  ,  C  =  (50  15  1  0) 

\-50 

-79 

-33 

-o) 

1 

\l/ 

A  model  of  order  nm  =  1  is 

Am  =  (  -0.576205 ) ,  i?m  =  (  1.073504 ) ,  Cm  =  (  0.588692 ) . 

This  model  yields  the  cost  J  =  0.104740.  A  model  of  order  rim  =  2  is 

_/ -0.532330  -0.598751\  „  _/  1.031824  \  0.588704 \ 

3.800771  -4.815122 /  ’  ~  V -3.103263  /  ’  ~  V  0.278923  J 

This  model  yields  the  cost  J  =  0.0269278.  A  model  of  order  rim  =  3  is 


/ -0.520312 

-0.731867 

-0.162146\ 

Am  =  2.888921 

-2.235622 

-3.721286  , 

V  -1.084500 

6.305395 

-0.746729  / 
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=  (1.020109  -2.114532  1.222072),  Cm  =  (0.586461  0.307967  0.105043), 

This  model  yields  the  cost  J  =  0.00148438. 

Example  9  [16].  The  system  is  given  by 


-6.2036 

15.054 

-9.8726 

-376.58 

251.32 

-162.24 

66.827 

0.53 

-2.0176 

1.4363 

0 

0 

0 

0 

16.846 

25.079 

-43.555 

0 

0 

0 

0 

377.4 

-89.449 

-162.83 

57.998 

-65.514 

68.579 

157.57 

0 

0 

0 

107.25 

-118.05 

0 

0 

0.36992 

-0.14445 

-0.26303 

-0.64719 

0.49947 

-0.21133 

0 

0 

0 

0 

0 

0 

376.99 

0 

'89.353 

376.99 

0 

0 

0 

0 

.  0 


0  ' 
0 
0 
0 
0 

0.21133 
0  > 


'0  0  0  0  0  1 
0  0  0  0  0  0 


1  0\ 

0  ij- 


A  model  of  order  rim  =  1  is 


Am  =  (-0.199272),  Bm  =  (0.631300  -0.00187918),  Cm  =  ( -0.187347  -354.430393), 

This  model  }rields  the  cost  J  =  27632.2.  A  model  of  order  rim  =  2  is 

_/ -0.199608  -0.0763006  \ 

“  V  3.331193  -13.275827  j  ’ 


„  _  /  0.631832  -0.00191612\  ^  (  -0-20105 

~  V -5.151821  -0.101952  /’  V -354.4141 

This  model  yields  the  cost  J  =  23262.3.  A  model  of  order  rtm  =  3  is 


_/ -0.199608  -0.0763006  \ 

“  V  3.331193  -13.275827 )  ’ 

00191612\  r  -  (  -0-201C 
1.101952  /’  V -354.414 


-0.201050  0.800899 

-354.414137  -66.187284, 


/  -0.198769  0.235666  -0.02481363 

Am  =  I  -1.087392  -0.912444  9.201811 

\  -0.115288  -9.502428  -0.0261157 


/ -0.630503  0.00216112 

Bm  =  (  -1.350879  -0.00377142 
V  -0.222387  -0.0526803 


_  (  0-2 
~  ^354. 


291338  -0.0265117  -4.035696 

1.222032  -164.479031  26.635498 , 


This  model  yields  the  cost  J  =  0.673079.  A  model  of  order  rim  =  4  is 


Am  — 


-0.198769  0.235667  -0.0248136  0.000915746 

-1.087390  -0.912440  9.201811  -0.00904508 

-0.115288  -9.502427  -0.0261155  0.00159031 

-5.465132  -11.698410  -1.929974  -37.554401 


5m  = 


-0.630503  0.00216112 

-1.350876  -0.00377141 
-0.222386  -0.0526803 
-8.666510  -0.0203036 


C*  — 


0.291340  354.222032 

-0.0265302  -164.479038 
-4.035692  26.635453 

0.0861885  -0.815898 
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This  model  yields  the  cost  J  =  3.22  •  10“^. 

For  this  example  with  rim  —  3,  4,  the  columns  of  the  initial  Jacobian  matrices  are  so  badly 
scaled  that  the  numerical  lineskr  algebra  in  HOMPACK  fails.  Modifying  HOMPACK  to  use  the 
LINPACK  subroutine  DQRDC  for  the  QR  factorization  of  the  initial  Jacobian  matrices  enables 
HOMPACK  to  successfully  overcome  the  ill  conditioning  and  find  a  solution. 

Table  1  gives  the  CPU  times  in  seconds  and  the  number  of  steps  needed  to  obtain  the  results 
for  each  example.  The  CPU  times  are  for  a  DEC  station  5000/200,  using  double  precision,  IEEE 
arithmetic,  and  the  MIPS  RISC  f77  compiler.  Table  2  gives  the  compzirison  of  the  optimal  projection 
equations  approach  and  the  input  normal  form  formulation  for  Examples  8  and  9.  The  asterisks  in 
Table  2  indicate  cases  that  required  special  numerical  linear  algebra  techniques  to  deal  with  severe 
scaling  errors. 

Table  1.  Algorithm  measures  for  input  normal  form  homotopy. 


example 

steps 

time  (sec) 

1 

1 

5 

0.06 

2 

1 

21 

0.13 

3 

1 

19 

0.10 

4 

1 

12 

0.14 

4 

2 

7 

0.20 

5 

1 

10 

0.12 

5 

2 

10 

0.22 

7 

1 

11 

0.22 

7 

2 

8 

0.30 

7 

3 

6 

0.46 

8 

1 

10 

0.20 

8 

2 

18 

0.50 

8 

3 

10 

0.65 

9 

1 

11 

0.71 

9 

2 

123 

8.0 

9 

3 

6 

1.3 

9 

4 

6 

1.9 

Table  2.  Comparison  of  methods 
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As  shown  by  Table  1,  the  input  normal  form  homotopies  can  be  very  efficient.  Also  there  is  no 
need  to  adjust  any  parameter  to  achieve  this  efficiency  (  although  to  obtain  the  minimum  solution 
of  Example  1,  some  adjustment  of  the  initial  point  was  necessary).  However,  note  that  the  potential 
ill  conditioning  of  the  input  normal  form  formulation  can  result  in  failure  (Example  6)  or  the  need 
for  extraordinarily  delicate  linear  algebra  (Example  9). 
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